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Abstract 

The ac Josephson effect in liybrid systems of a normal mesoscopic conductor 
coupled to two superconducting (S) leads is investigated theoretically. A gen- 
eral formula of the ac components of time-dependent current is derived which 
is valid for arbitrary interactions in the normal region. We apply this formula 
to analyze a S-normal-S system where the normal region is a noninteracting 
single level quantum dot. We report the physical behavior of time-averaged 
nonequilibrium distribution of electrons in the quantum dot, the formation 
of Andreev bound states, and ac components of the time-dependent current. 
The distribution is found to exhibit a population inversion; and all Andreev 
bound states between the superconducting gap A carry the same amount of 
current and in the same flow direction. The ac components of time-dependent 
current show strong oscillatory behavior in marked contrast to the subhar- 

monic gap structure of the average current. 
74.50.+r, 73.40.Gk, 73.23.-b, 72.15.Nj 
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I. INTRODUCTION 



Quantum transport properties of mesoscopic conductors coupled to two superconducting 
(S) leads have been extensively investigated in the last decade both theoretically and exper- 
imentally The mesoscopic conductor in question is usually not a superconductor itself, 
but it can be a quantum point contact (QPC) a quantum dot (QD) a tunnel 

barrier, a normal metal |TO|JTl[l , and even a molecule such as a nanotube [p!2| -[T^. The physics 
of these hybrid device structures, in the form of S-normal-S, has profound implications to 
both fundamental understanding of quantum transport at reduced dimensionality and to 
practical applications in nanoelectronics. 

One of the main transport characteristics of a S-normal-S device structure is that particles 
in the normal region can undergo multiple Andreev reflections by the two superconducting 
leads. If the normal region is ballistic, a consequence of the coherent superposition of these 
multiple Andreev reflections is the formation of Andreev bound states [|l|,|l^. The Andreev 
bound states are important because they carry current including the supercurrent. On the 
other hand, if the normal region is diffusive, a so-called supercurrent-carrying density of 
states, instead of the Andreev bound states, gives the ability for carrying supercurrent [jl6l . 



The multiple Andreev reflection is also known to generate subharmonic gap structure in the 
behavior of Jq = /o(^)) where Iq is the average current and V is the bias voltage P-p|,p!0|. 
More recently, the subharmonic gap structure is used to measure transmission probability 
of each channel in a multi-channel QPC device |T^,[TB[ . 



Another important and interesting transport characteristic of S-normal-S devices is the 
Josephson effect which gives rise to a dc supercurrent at zero bias, and an ac current at non- 
zero bias. Previous theoretical analysis have focused on the dc Josephson effect at zero bias 
0, and the subharmonic gap structure of the average current at a non-zero bias |]3|-|8|,p^ 



However, the ac Josephson effect, which arises at a non-zero bias, produces a current that 
is a function of time t. Therefore it is an important task to theoretically understand the 
time dependent current in addition to understanding its time-average. To the best of our 
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knowledge, so far there have been only two works which involve a time-dependent current 
for S-normal-S devices. Cuevas et. al. have investigated the ac component of the 
time-dependent current for a S-QPC-S system [|]. Bratus et. al. have investigated the 
time-dependent current in a S-quantum-constriction-S system by considering an arbitrary 
normal electron transparency and discussing the current property at the small bias limit. 
In both these works, the normal region of the S-normal-S device was simplified so as not to 
have any electronic structure: it is simply described by a constant transmission coefficient 
which is independent of energy e. Given the interesting physics already discovered by these 
previous investigations, it is indeed not difficult to expect that even richer physics would 
arise if the normal region has its own electronic structure. 

It is the purpose of this work to further investigate the ac Josephson effect in S-normal-S 
device systems, and we focus on issues not resolved by the previous analysis. In particular, we 
consider a mesoscopic S-normal-S device with an arbitrary normal region which may have 
its own electronic structure and/or strong electron-electron interactions: for this general 
situation we have derived the expression of the ac current. As an application we then 
investigate a specific case in which the normal region is a ballistic quantum dot having 
a noninteracting single energy level, for which we investigate the intradot nonequilibrium 
distribution of electrons, the local density of state (LDOS), and the time-dependent current. 
Our main findings are: (i) The intradot electronic distribution shows a population inversion 
property. This property is distinctly and qualitatively different from that of the case where 
the normal region is diffusive, (ii) At small bias voltages such that eV < A where A is the 
superconducting gap, a series of Andreev quasi-bound states is found to emerge within the 
gap. Their weights are not the same but they carry equal amount of current in the same 
direction, as well, their electronic occupations are all 1/2. This is qualitatively different 
from that of the zero bias case in which the successive Andreev bound states carry opposite 
current, (iii) The ac current component versus bias V shows an oscillatory behavior. The 
amplitude of oscillation of the nth component is largest at about V = A/n. At small bias, 
the high-order components quickly increase, and the time-dependent current versus time t 
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deviates from a sine-like curve. 

The rest of the paper is organized as follows. In Sec. II, the model Hamiltonian is 
presented and a general formula for the ac current component is derived. In Sec. Ill, 
ac Josephson effect for a simple S-normal-S device with a noninteracting normal region is 
investigated. The intradot electronic distributions, the Andreev quasi-bound states, and ac 
current components are presented in this Section. Finally, a brief summary is given in Sec. 
IV. 



II. MODEL AND FORMULATION 



We assume the S-normal-S device system to be described by the following Hamiltonian: 



where 



H = ^ Ha + Hcen + Ht , (l) 
a=L,R 



k,a k 



k,j,a,a 

Ha (a = L^R) describes the left/right BCS superconducting lead with the superconducting 
energy gap A^. H^en is the Hamiltonian of the normal region of the device, and c]^(cjo-) 
are the creation (annihilation) operators of an electron in state ja of the normal region. 
Hint models interactions in the normal region whose form depends on specific physics prob- 
lems under consideration. In this Section we consider the general case without specify- 
ing its concrete form. In deriving the formula for the transport current, we permit the 
device normal (central) region to have various interactions, such as the electron-electron 
Coulomb interaction, Uj(Tvh^i^]<7^i<^^]icyx'^hoi'^ electron-phonon interaction, 

MjgCj^Cja{d^q + d_q) + '}2T^^qd\dq] the tuuneliug coupling between different states of the 



normal region, J2 tijcl^Cj^ + H.c. ; and so on. Ht of Eq.(|T]) denotes the tunneling 
Hamiltonian between the superconducting leads and the normal region of the device, and 
taj is the hopping matrix. In order to obtain the Hamiltonian (1), we have performed a 
unitary transformation, then the superconducting initial phase (pa and the terminal voltage 
Va emerge in the Hamiltonian Ht- |T9| , |20[] 

The total current of superconducting lead a {e.g. a = L) flowing into the device normal 
region can be calculated from evolution of the total number operator of electrons in that 
lead, A''^ = J2k,a ^Lka'^LkcT- Then we have (in units of h = 1): PU|-P^ 

hit) = -e< iVi(t) >= te < [Nl, H] > 



2eRe ^ tue 

k.i 



(5) 



where 



Gfrkit,tl) = ^ 



is the distribution Green's function in the 2x2 Nambu representation, and is the Pauli 
matrix. In this paper, we use the notation that "A" means quantity A to be a 2 x 2 matrix. 
To proceed we need to solve the Green's function Gfj^j^it^t). We assume that the leads 



do not have any interactions except the quadratic pair potential correlation, we have: P0|j22 



(6) 



where g^j^{ti,t) is the exact Green's function of the left superconducting lead, [|,|T^] and 
^Ljiti) is a 2 X 2 hopping matrix defined by: 



iUt) 



tLjC V 



■eV,t 







-the 



(7) 



and Gfjit,ti) are the retarded and distribution Green's functions in the device 
normal region. They are defined by: 
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GUt,h) = -i9{t-h) 



^ < {Ci^{t),c]^{ti)} > < {Ci^{t),Cji{ti)}>^ 



(8) 



Gf, {t, ti 



(9) 



Substituting C^j^f^it.t) into Eq.(5), assuming tij is real, the current /L(t) can be expressed 
in terms of the Green's functions of the device normal region, as: 



h{t) 



-2elm 



J 27r 



(10) 



where / 



l/r[ 



is the Fermi distribution function of electrons in the left/right 



superconducting lead. is defined as: 



for Al < |e|. PL(e) = i?e[/?z.(e)] 



/?L(e) 
|e|-A)- 



for > |e|, and /^^(e) 



is the dimensionless BCS 
density of states, i.e. the ratio of the superconducting density of states pf (e) to the normal 



density of states (e). F is the linewidth matrix function defined by Ti-ij = 2ntiit*^jp^ (e 



in which we have assumed that F^ is independent energy e. |2J] In this paper, we use boldface 
letters to denote quantities representing matrices whose matrix elements are calculated using 
states i,j of the device normal region. Finally, is a compact notation. 



/ 



SL(e) 



-ieVLih-t) 



_^^-i4>L-ieVL{ti+t) 



V 



(11) 



_^^i't>L+ieVL(ti+t) ^ieVL{ti-t) 
e 

The formula Eq.(lO) describes the current using Green's functions of the normal region. It is 
a general formula and can therefore be applied to situations involving arbitrary interactions 
in the normal region and is also applicable at nonequilibrium {e.g. at a high bias V^). If the 
normal region is coupled to multiple superconducting leads or to some extra normal leads, 
Eq.(lO) is still vahd. 



In the following we fix Vl = |^ so that the left superconducting lead is taken as the 
potential ground, then T,^ reduces to: 



(l 



— e 



(12) 



Note that the superconducting phase difference between the two leads is a time dependent 
periodic function with a period T = 271 /u, where u = 2eV and V = Vl — Vr is the bias 
voltage between the leads. Therefore the time- dependent current /l(^) is also a periodic 
function with the same period T because the Green's functions have the property G{t, ti) = 
G{t + T,ti + T). |^6|] Then we can take the conventional Fourier expansion for the current 
hit): 



inuit 



lLit)=J2lLne 

n 

and take the double Fourier expansion for the Green's function: 

de 



(13) 



G{t,t,) = ^e*«-*i| _le-(*-*i)G'„(e 



(14) 



To simplify notation in the following analysis, we introduce quantities Gmn(e) = Gn~m{^ + 
muj) and Xi(t), 

de 



X^{t) = -2e / dti I —e'<^-^^^Tr 



Pi(e)^(e)G^(t,ti)+/32(e)G<(t,ti) TlJ^l&A , (15) 



so that hit) = Im[lL{t)]. 

Then the Fourier component of ac current is obtained as: 



1 



Ln) 



(16) 



and 



Ln 



I de ^ 
-2e / —Tr 

2lT 



■nO\ 



o/?2(e)G: 



■nO\ 



(17) 



Eqs. (0,113) are the first central results of this work. They describe ac components of the 
time-dependent current of a S-normal-S device system in terms of the Fourier component of 
the Green's function G!l„o(e) and Gf„o(^) of ^^e normal region. These formula, Eqs. (13), 
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(16), and (17), are valid for arbitrary interactions the normal region may have, for nonequi- 
librium situations, and for devices with other normal leads. They can not, however, be 
applied to devices with more than two superconducting leads. 

When bias voltage V is zero the current Iiit) is independent to time t, then the current 
reduces as: 



-2elm 



de 
2^' 



-Tr 



^(e)p,(e)G^(e) + -/52(e)G<(e) 



(18) 



III. NONINTERACTING NORMAL REGION 

In this section we apply the general expressions for the ac current derived above to an ex- 
ample of a S-normal-S device where the normal region has no electron-electron interactions. 
For this situation, the Hamiltonian i^cen can be written as: 

«J,o-(i>i) 

= Y.Hcen,a • (19) 

This Hamiltonian describes a multi-level noninteracting quantum dot for which tij = 0. It 
also can describe a typical tight-binding lattice model, in which tij ^ 0, the second term in 
Eq.(p!9D denotes the coupling between different lattice sites. 

For the specific Hcen of Eq. (P^ , we can solve the Green's functions G^„(e) and G^„(e) 
using the Dyson equation and the Keldysh equation: G'' = g'' + G^S'^g'^, and G^ = 
G^'S^G'^. Here g*" is the exact Green's function for the device normal region without 
coupling to the leads, and it can be easily derived as: 



f) piHcETilit-ti) 

V / 



(20) 



S'' and S"^ are the retarded and distribution self-energies due to coupling to the leads, with 
S^«)(t,ti) = E?<)(t,ti) + ti<\t,t,) and 



^ / ^^L{R),ijfL{R){e)pL(R){^)e ''^h{R) . 



(21) 



(22) 



The Fourier space form of these quantities are easily obtained (notice that Vl — and 



3mn 



^mnl (Cm " '^cerA + «0+) 



^L;mn(^) " " 2 r'L<^mn/3L(er„)EL(e^ 



^mnPRi^m+i) 



5^,n+i/5i^(e^_i)f^e^<^« <^mn/3il(e^ 1) 



^L;mn(e) = L^mnf L{em)pL{em)^L{er 



5mn/L(e^+i)pil(e^+i) 5m,n-lfL{em+l)pR{^m+\)7^(^ "''"'^ 



^m,n+i/L(e^_i)pR(e„^_i)7^e'*« 5mn/i?(em-i)PR(em-i! 



(23) 
(24) 
(25) 

(26) 
(27) 



where — e-\-xuj. Similarly, the Fourier space form of the Keldysh equation and the Dyson 
equation are: 



h,l2 



(28) 
(29) 



If G'[^^{e) has been solved, then from the Keldysh equation (28), G^„(e) can be obtained 
straightforwardly. Therefore in the following we only need to solve the retarded Green's 
function G^„(e). 

From the Dyson equation (29) we have: 



^mn;ll ~ &mn;ll"mn "T '-^mn;ll^nn;116nn;ll '^m/;12^in;21&nn;ll ' 

I 

^mn;12 " '-^mn;12^nn;22&nn;22 ^mi;ll^in;126nn;22 ' 



(30) 
(31) 



where we have suppressed the argument e. From Eq.(31), one has: 



^mT?,:12 ~ y^^mi!:ll^[w:12 r--1 _ -^r " ("^^) 

I &nn;22 ^nn;22 



Substituting this expression to Eq.(30) one easily finds: 



Gr _ ^rnn i_ \^ r<r i ) /qqN 
mn;ll " r-1 _ vr ^ 2^ ^ml;ll^ln : W-^J 

&nn;22 ^nn;22 I 



where 



1 1 

Bmnl^) = y^,'^ml:12~T^ ^ ^'ln;2l'~F^ ^ ■ (34) 

I Sll;22 ^ll;22 6nn;ll ^nn;U 

Note Bmn 7^ only when m = n, n ± 1. The quantity has a clear physical meaning: 
it describes the intensity of Andreev refiection processes, an example is shown in Fig.l in 
which a particle in the normal region undergoes twice Andreev reflections. Then by iterating 

Eq.(33), G'^^ .^j^ can be formally solved, 

^mn;ll ~ r-1 _ yir r-1 _ yr ^mn , (35) 

6nn;ll nn;ll omm;ll mm;ll 



where 



= B^„ + ^ B^/ Y;„ . (36) 
I 

Similarly, the quantity Y^„(e) has a clear physical meaning: it gives the intensity of the 
process for which an electron having initial energy e + ncu ends up with flnal energy e + imu 
after going through all possible multiple Andreev reflections in the normal region. Eq.(36) 
can only be solved numerically and after Y^n is solved, from Eqs.(35) and (32) G^^.^j^ and 
^mn;i2 ^au be obtained immediately. Finally, Gmn;2i ^mn;22 can also be calculated 
using following equations which are derived from the Dyson equation: 

Gfmn;21 = X] TF^l _ -^r ^ ml;21^ ln;ll ' (37) 

I 6mm;22 ^mm;22 

r -1 
f^r I \ ^ yir f-^r 

^mn;22 — r-1 _ yr „r-\ _ yr ^ ml\2V^ ln\\2 ■ y^°) 

6mm;22 '^mm;22 ; &mm;22 ^mm;22 
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With GJ„„ and solved, from Eq.(17) the ac component and time-dependent current 
can be calculated without further complications. 

In the rest of this section, we present numerical results for which some further simplifica- 
tions are made. We reduce the device normal region to a quantum dot with a spin degenerate 
single level, i.e. Hcen = Y^^dC^a^a- For this case the boldface matrices reduce to a C num- 
ber. We also take A = A/, = A^j = 1 as the energy unit and only consider devices with 
symmetric barriers {Tl = Tj?)- It should be mentioned that since we have assumed a spin 
independent intradot level and hopping elements tL{_R), < cj(ti)c|(t) > should be equal to 
< c{(ti)q(t) >. Following this we have Gf^{t,ti) + ^^2(^1^^) = - [^11(^.^1) - 
and G<(t,ti) = -[G<(ti,t)]^ The Fourier forms are a<„.n(e) + G<^^_^.^^{-e) = 
— [G'J„„.]^;^(e) — GJ^m-iil^)] ^'^d G^^{t) = — [G'J^„(e)]"''. These relationships provide very strong 
checks on our analytical derivations and numerical calculations which we present in the 
following subsections. 

A. Intradot distribution of electrons 

In this subsection we present results of the intradot distribution of electrons for the 
S-normal-S device. Because of the finite bias voltage V, the current, intradot occupation 
number of electrons, local density of states (LDOS), and the intradot distribution of elec- 
trons, are all functions of time t. The time average occupation number of electrons on the 
intradot state t is: (same for state |) 

^G'<o;n(e) • (39) 

The integrand of (|3|), =^G< 

■11(f)) is the time-averaged occupation number of electrons with 
energy e. Here, subscript "11" are indexes of the 2x2 Nambu matrix element, and "00" are 
indexes of Fourier component. The average LDOS is given by LDOS{e) = —^Im[GQQ.ii{e) + 
^oo-22(~f)]- The average intradot distribution of electrons can be obtained from the average 
occupation number at energy e and the average LDOS{e) p^ . 
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It is important to emphasize that the distribution of electrons can be experimentally mea- 
sured [p9| , |30| . For example, recently Pierre et. al. have measured this distribution for 



a S-normal-S device where the normal region is a diffusive mesoscopic metallic wire. They 
reported a multiple step structure for the distribution of electrons in that device |30(| . 

Fig. 2 shows the average intradot distribution of electrons at different bias voltage V for 
our system with a very large coupling F. When F is large, coupling between the supercon- 
ducting leads and the normal region is strong, therefore the device behaves like a S-ballistic- 
normal-conductor-S system. The property of the electron distribution in this situation is 
the following. When miniVL — A,Vr — A) < e < maxiVL + A, Vr + A), the distribution is 
a constant, i.e. fd{^) ~ 1/2 for symmetric couplings. When e goes away from this region, 
the distribution quickly rises (or drops) to unity (or to zero) for e < min{VL — A, Vr — A) 
(or for e > max{VL + A, Vr + A )). 

To contrast with the experimental results of Pierre et. al. |^, here the distribution is 



a constant instead of the multiple step structure between the gap, even though multiple 
Andreev reflections do occur in our system. This difference originates from the different 
property of the central device region, i.e., our normal region is ballistic while that in Pierre 
et. al. experiment is diffusive ^0|. In order to explain it in more detail, the inset of Fig. 2 
shows a particular multiple (two) Andreev reflection process. To start, an incident electron 
at Ai below the gap of the left lead tunnels into the QD, it passes two Andreev reflections 
(through the points labelled as A1-A6) inside the QD and finally tunnels into the right 
lead (at A^) which is higher than the gap of the right lead. Due to the ballistic nature of 
the QD, the distribution of electrons at point Al is the same as at A2, the distribution of 
holes at A3 is the same as at A4, while distribution of electrons at A5 is the same as at 
A6. When F is large, the probability of Andreev reflection inside the QD within the energy 
gap is very close to unity |^, and hence the hole distribution at A3 is, to a very good 
extent, the same as the distribution of electrons at A2. Similarly the hole distribution at 

12 



A4 is approximately the same as the electron distribution at A5. We hence conclude that 
for the ballistic normal region, the distribution of particles (electrons and holes) along this 
path is the same everywhere, except at the abrupt change during the tunneling process at 
Ai and from and to the two leads. Moreover, for symmetric barriers, the distribution 
of particles along the A1-A6 path must be 1/2. This explains why we obtained a constant 
1/2 distribution at mm(Vi — A,Vr — A) < e < maxiVL + A,Vr + A) as shown in Fig.2. 
This also explains why we expect a different distribution when the normal region is diffusive: 
for a diffusive conductor the distribution at Al and A2 must be different due to diffusive 
scattering between the two points, therefore the distribution of particles will continuously 
vary from one to zero along the path A1-A6. 

Next, we investigate the distribution of electrons for F ~ A, the results are shown in 
Fig. 3. For this case, a most prominent behavior of /d(e) is that it oscillates as a function of 
e. The oscillations also become more rapid when bias voltage V is reduced. An oscillatory 
/d(e) means its value is not necessarily smaller for larger e, hence a "population inversion" 
is possible. This population inversion originates from the non-monotonic probability of 
Andreev reflections. For example, /d(e) has a dip at e = — A, due to the following 
reason. For an incident electron coming from the left lead with energy Vr — A, this electron 
has a small but nonzero probability to pass the left barrier. After tunneling through, it 
reaches the right barrier where an Andreev reflection occurs. Because this electron has 
energy e = — A, the Andreev reflection occurs with probability one |3l|. Therefore the 



distribution of electrons at this energy e is very small. When e deviates from Vr — A, the 
probability of Andreev reflection decreases leading to a larger f^, hence we expect a dip in 
fd to emerge at e = Vr — A. 

B. Local density of states 

In this subsection, we investigate another important quantity, the LDOS. We will mainly 
discusses Andreev bound states at a finite bias V. If bias V > 2A, multiple Andreev 
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reflections are very weak hence no Andreev bound states can form in tlie QD. In tliis case 
tlie intradot level is only slightly shifted due to a non-zero real part of the self-energy S^, 
the level half-width is still on the scale of Tljh-, and an extra structures (a dip and a peak) 
emerge in the curves of LDOS{e) versus e at the superconducting gap (not shown in here). 

Much more interesting is the case of V < A, shown in Fig. 4 at different bias V . A series 
of very narrow peaks emerge in LDOS{e), clearly indicate the formation of Andreev bound 
states inside the QD. Note that they are not rigorous bound states but are quasi-bound 
states with a finite life time, and after many Andreev reflections the particle can leave the 



QD. This is different from the zero bias situation |jT5[. The half-width of Andreev bound 
states is much narrower than F. With a decreasing bias V, they become even narrower with 
a higher intensity. The average distance between two successive Andreev bound states is 
about eV. When neV and (n + l)eV {n = 0, ±1,±2,...) are within the gap, there exists 
an Andreev bound state between e = neV and (n + l)eV. Moreover, these Andreev bound 
states are symmetrically distributed at the two side of Vl and Vr. This means the following: 
when an incident electron from below the gap aligns perfectly with an Andreev bound state 
of the QD, even after many Andreev reflections it will always stay on the Andreev bound 
state until it leaves the QD (see inset of Fig.4(a)). Along this path, the particle goes through 
all Andreev bound states, and a resonance multiple Andreev reflection occurs. Occasionally, 
a quenching of Andreev bound state is observed to occur. In this case, a specific Andreev 
bound state may have very low LDOS at a specific bias V, an example is indicated by the 
arrow in Fig.4(b). 

The results of Fig. 4 is obtained by fixing the intradot level to zero {i.e. at the center 
of the gap). Next, we investigate how are Andreev bound states affected when 7^ 0, the 
results shown in Fig. 5. With 7^ 0, the Andreev bound states are shifted in their positions, 
but their physical characteristics are the same as those of = 0. The amount of shift is 
not ed but much smaller and two successive Andreev bound states are shifted in opposite 
directions. If an Andreev bound state is in the energy range from e = neV to {n + l)eV, it 
stays in this range at any value of e^. Their heights vary with e^, when is in the range of 
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neV to {n + l)eV, the peak in this range reaches a maximum value. 

An important property of the Andreev bound states is their abihty to carry current. 
From Eqs.(16) and (17), the time-averaged current density jo(e) is obtained to be: 



Me) 



-^ImTr I [h{e)pLie)Gl,{e) + ^/3I(6)G<o(e)] T^lct.] . (41) 



The current density is related to time-averaged current as Jq = / (iejo(e)- In Fig.6, we 
show intradot distribution of electrons fd (solid curve in Fig.6a), LDOS (dotted curve in 
Fig. 6a), and the time-averaged current density (Fig.6b) jo(e)- Several observations are in 
order, (i). Although /^(e) is oscillating between and 1 in a complicated manner, its value 
at each Andreev bound state (the peak positions of the dotted curve) is alway 1/2. This 
is because resonant multiple Andreev reflections occur along the path of Andreev bound 
states (as shown in the inset of Fig. 4(a)). (ii). The current density jo{e) is dominated by 
a series of peaks located precisely at the energies of Andreev bound states. This is a clear 
indication that current is carried by Andreev bound states. When mm(VL — A,Vr — A) < 
e < maxiVi + A, Vr + A), the peaks of jo{e) all have the same height: this means each 
Andreev bound state carries exactly the same amount of current in the same flow direction. 
The reason for this peculiar behavior is simple. Along the path of Andreev bound states 
(inset of Fig.4(a)), all the electrons move in one direction while all the holes move in opposite 
direction, and along any one path the particle current must be same everywhere. Therefore 
the Andreev bound states carry same amount of current in the same direction. This property 
is qualitatively different from that of the zero bias case , in which the successive Andreev 
bound states carry current with opposite sign. 



C. The current 

The time-averaged current Iq of S-normal-S systems has been extensively investigated 
both theoretically and experimentally. A main characteristic of the I-V curve Io{V) is its 
subharmonic gap structure at = 2A/n our results are shown in Fig. 7. The I-V 
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curves also exhibit subharmonic gap structure with a concomitant appearance of negative 
differential conductance. These results are in agreement with those reported recently by 
Yeyati et. al. and Johansson et. al. In the following, we focus on the ac component 



of the current. 

From Eqs.(13) and (16), we decompose the time-dependent current into its dissipative 
contribution J^, and nondissipative contribution 

h (t) = h + Yl ^Lo cos nujt + IIq sin nut , (42) 

n n 

where /£„ = Im(lLn +lL-n) and /£„ = ReiTin — Il-u)- Fig. 8 and Fig.9 show the first and 
second ac components of and /£„ as a function of bias V , and they are marked by a strong 
oscillatory behavior. The period of oscillations is roughly given by which is dependent 
on bias V. Generally, for < eV < (m=l,2,...), we found that the ac components 
oscillate from a maximum to a minimum or vise versa. When > the components 
/£„ and /£„ quickly decay to zero. When eV ~ ^, the amplitudes of the oscillations reach 
maximum. At eV 0, J£„ decays to zero while J£„ keeps a finite value. These behaviors are 
different from those devices whose normal region has no electronic structure. For instance, 
the result of S-QPC-S system shows no oscillation [Q. 

The time- dependent current /^(t) is shown in Fig. 10. /l(^) is a well known oscillatory 
function of time t with a frequency uj = 2eV. When bias V is large, eV > A, the high- 
order Fourier components have negligible contribution and /l(^) can be approximated by 
Iiit) ~ /q + lLisin{ijjt + (f)). On the other hand, when V is small, high-order components 
contribution substantially and /l(^) deviates from a simple sine-like curve. 



IV. CONCLUSIONS 

In this work, we have derived a general formula for ac components of the time-dependent 
current of arbitrary ballistic S-normal-S systems where the normal region has its own elec- 
tronic structure. The formula (Eq.(17)) is valid even when there is a strong interaction in 
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the normal region of the hybrid device. We then apphed this result to study ac Josephson 
current for a system with the normal region being a noninteracting single level quantum dot. 
The average intradot distribution of electrons, the average intradot density of states, and 
ac components of the time-dependent current are investigated in detail. The distribution 
exhibits an interesting population inversion, a result that is qualitatively different from that 
of the diffusive normal region. A series of Andreev bound states are formed at bias V < A in 
our system. The peak heights of LDOS for these Andreev bound states are not the same, but 
each state carries the same amount of current. The distribution of electrons at the Andreev 
bound states are all the same, e.g. equals 1/2 for symmetric tunnel barriers. In general, 
the ac components of the time-dependent current has an oscillatory behavior against bias. 
Depending on the value of bias, the high-order ac components may or may not contribute 
to the total time-dependent current, leading to a non-sine- like or a sine-like dependence on 
time for the total current. 

Finally, we comment on the eV limit for the S-QD-S system of this work. While our 
general current formula, Eq.(17), is valid for any bias, how to correctly include important 
physical factors in an actual computation of the various quantities of Eq.(17), needs to be 
discussed. When bias is very small, eV <^ A, an incident electron from below the gap of 
the left superconducting lead undergoes many Andreev reflections in the QD so as to go 
above the gap of the right superconducting lead before exiting the QD. Therefore the dwell 
time Tp of the particle in the QD becomes long. At the limit eV —>■ 0, Tp tends to large 
values. When Tp is larger than the mean inelastic scattering time, the intradot relaxation 
effect should be considered in calculating the Green's functions involved in Eq.(17). When 
there is no electronic structure in the normal region of the device, for instance in a S-QPC-S 
system @J^, the eV limit has a variety of different regimes depending on an inelastic 
scattering rate parameter 6 and a transmission probability of the QPC 0,0. For our S- 
QD-S system, while relaxation in the leads can similarly be included by introducing the 
same parameter 6 into the Green's function of the leads 0], this simple phenomenological 
approach can not be applied in the normal QD region. This is because distribution of 
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leads is determined by their chemical potential, however the distribution in the QD must 
be calculated self-consistently for our system. Indeed, if one introduces a finite 5 in the QD 
Green's function, current conservation will be violated. A proper treatment of this problem 
is, perhaps, to explicitly introduce an electron-phonon interaction term in the Hamiltonian. 
This is a very complicated problem to solve and we hope to be able to report such an analysis 
in the future. 
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FIGURE CAPTIONS 



Fig. 1 A schematic diagram for the transport process consisting of two Andreev reflections, 
(a). The particle is first Andreev reflected by the left superconducting lead, then it is 
by the right superconducting lead. This is described by the quantity Boi(e). After this 
process, the particle energy reduces by 2eV {i.e. u — 2eV). (b). The particle is first 
Andreev reflected by the left (right) lead, followed by another reflection at the same 
lead. This process is described by quantity Boo(e). After this process, the particle 
energy does not change, (c). The particle is first Andreev reflected by the right lead, 
then by the left lead. It is described by quantity Bo,_i(e). After this process, the 
particle energy rises 2eV. All processes with an even number Andreev reflections can 
be decomposed to the three processes plotted here. All processes with an odd number 
of Andreev reflections can be decomposed to the even case plus one more reflection. 

Fig. 2 The time-averaged intradot distribution of electrons versus energy e at large F, = 
Fr = lOOOA. Temperature KbT = 0.05A, = 0, d = {6 is the inelastic scattering 
rate in two superconducting leads), and 4>l = 4>r = 0- Note the fact that the time- 
averaged distribution, LDOS, and the ac components of the current are all independent 
with initial values of (j)^ and (p^at 5 — 0. Inset: schematic diagram showing a multiple 
(two) Andreev reflection process. 

Fig. 3 The time-averaged intradot distribution of electrons versus energy e at general QD 
parameters, F^, = F^^ = 1.5A. Other parameters are the same as those of Fig. 2. 

Fig. 4 The time-averaged LDOS versus e at different bias V. KbT = 0.1 A, Tl^Tr^ 0.8A, 
ed — 0, and 5 — 0. The downward arrow in (b) points to an Andreev bound state 
with a very small LDOS. Inset in (a): schematic diagram showing a multiple Andreev 
reflection which passes through the Andreev bound states indicated by the thick solid 
lines in the QD. 
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Fig. 5 The time-averaged LDOS versus e at different level positions e^. Vr — — 0.3A and other 
parameters are same as those of Fig.4. Different curves correspond to = 0.15 A, OA, 
— 0.15A, — 0.30A, and —0.45 A, along the arrow direction. 

Fig. 6 (a). The time-averaged LDOS (dotted) and the time-averaged distribution of elec- 
trons (solid) versus e; (b) the time-averaged current density, ed = —0.15 and other 
parameters are same as those in Fig. 5. 

Fig. 7 The time-averaged current Iq versus bias V at different F. Other parameters: KbT — 
O.IA, ed = 0, 5 = 0.005A, 0^ = 0ii = 0. 

Fig. 8 The dissipative ac components and versus bias V at different F. Other param- 
eters are the same as those of Fig. 7. 

Fig. 9 The nondissipative ac components and versus bias V at different F. Other 
parameters are the same as those of Fig. 7. 

Fig. 10 Time-dependent current /l(^) versus time t at different bias V. Tl = = 0.8 A and 
other parameters are the same as those of Fig. 7. The curves labelled 1 to 5 correspond 
to y = -Vr = 0.2A, 0.5A, l.OA, 1.5A, and 3.0A, respectively. 
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